
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# DISCLAIMER AND GENERAL INFORMATION
#
# File: App_A_main_results.R
# Purpose: Appendix file for main results
# Date: 10 July 2024
# Data: Survey data (pulled through 0_cbam_prep.R)
#
# Technical disclaimer:
# All analyses in R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"
# R Studio 2024.04.2 Build 764 ("Chocolate Cosmos" Release (e4392fc9, 2024-06-05) for Windows)
# Windows 10 Enterprise, 64-bit
# 12th Gen Intel(R) Core(TM) i7-1255U 1.70 GHz with 16GB RAM
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (A) Load data and packages ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Run source file to load data
source("0_cbam_prep.R")

# Load packages
library("tidyverse")
library("janitor")
library("ggpubr")
library("modelsummary")
options(modelsummary_format_numeric_latex = "plain")
library("marginaleffects")
library("xtable")
library("tinytable")

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (B) Regression table of main results ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: collapsed treatment conditions
ols_de <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany",]) 
ols_hu <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary",]) 
ols_ch <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland",]) 
ols_uk <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk",]) 

order_full <- c("factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)4"="Treatment: Both")  

models.full <- list(
  "Germany" = ols_de,
  "Hungary" = ols_hu,
  "Switzerland" = ols_ch,
  "UK" = ols_uk)


# Produce results table (Table A.1)
tab <- modelsummary(models.full,
             statistic = 'conf.int',
             conf_level = .95,
             coef_rename = order_full,
             gof_map=c("nobs"),
             title=c("Main regression results for CBAM support"),
             output="tinytable") %>%
             theme_tt("placement", latex_float="h!")

tinytable::save_tt(tab, "./main_cbam_support.tex", overwrite=TRUE)




# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (C) Main results for all policy frames ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: all conditions
ols_de <- lm(cbam_support~factor(exp_group), data=dta[dta$country=="germany",]) 
ols_hu <- lm(cbam_support~factor(exp_group), data=dta[dta$country=="hungary",]) 
ols_ch <- lm(cbam_support~factor(exp_group), data=dta[dta$country=="switzerland",]) 
ols_uk <- lm(cbam_support~factor(exp_group), data=dta[dta$country=="uk",]) 


order_full <- c("factor(exp_group)12"="Climate frame:\n Both",
                "factor(exp_group)11"="Climate frame:\n Prices",
                "factor(exp_group)10"="Climate frame:\n Jobs",
                "factor(exp_group)9"="Climate frame:\n Control",
                "factor(exp_group)8"="Trade frame:\n Both",
                "factor(exp_group)7"="Trade frame:\n Prices", 
                "factor(exp_group)6"="Trade frame:\n Jobs",
                "factor(exp_group)5"="Trade frame:\n Control",
                "factor(exp_group)4"="Both",
                "factor(exp_group)3"="Prices",
                "factor(exp_group)2"="Jobs")

order_full2 <- c("factor(exp_group)12"="Climate frame:\n Both",
                "factor(exp_group)8"="Trade frame:\n Both",
                "factor(exp_group)4"="Both",
                "factor(exp_group)11"="Climate frame:\n Prices",                
                "factor(exp_group)7"="Trade frame:\n Prices",                 
                "factor(exp_group)3"="Prices",
                "factor(exp_group)10"="Climate frame:\n Jobs",
                "factor(exp_group)6"="Trade frame:\n Jobs",
                "factor(exp_group)2"="Jobs")


models.full <- list(
  "Germany" = ols_de,
  "Hungary" = ols_hu,
  "Switzerland" = ols_ch,
  "UK" = ols_uk)

# Plot coefficients (Figure A.1)

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(models.full, coef_map=order_full2, coef_omit="Interc", conf_level = 0.95) + facet_grid(~model) +
  labs(x="Coefficients",
       y="Experimental conditions",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic() +
  scale_color_manual(values = cols) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none")
  #coord_flip()
  #theme(axis.text.x = element_text(angle = 90))
p1


# Export plot to PDF
ggexport(p1, filename="./Exp_support_byframe.pdf", width=10, height=5)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (D) Absence of framing effects ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Germany
de <- data.frame(Germany = 
                 # Framing effects for jobs treatment
                 c(hypotheses(ols_de, "b2=b6")$p.value[1], 
                 hypotheses(ols_de, "b2=b10")$p.value[1], 
                 hypotheses(ols_de, "b6=b10")$p.value[1],
                 # Framing effects for prices treatment
                 hypotheses(ols_de, "b3=b7")$p.value[1], 
                 hypotheses(ols_de, "b3=b11")$p.value[1], 
                 hypotheses(ols_de, "b7=b11")$p.value[1],
                 # Framing effects for compound treatment
                 hypotheses(ols_de, "b4=b8")$p.value[1], 
                 hypotheses(ols_de, "b4=b12") $p.value[1],
                 hypotheses(ols_de, "b8=b12")$p.value[1]))

# Hungary
hu <- data.frame(Hungary = 
                   # Framing effects for jobs treatment
                   c(hypotheses(ols_hu, "b2=b6")$p.value[1], 
                     hypotheses(ols_hu, "b2=b10")$p.value[1], 
                     hypotheses(ols_hu, "b6=b10")$p.value[1],
                     # Framing effects for prices treatment
                     hypotheses(ols_hu, "b3=b7")$p.value[1], 
                     hypotheses(ols_hu, "b3=b11")$p.value[1], 
                     hypotheses(ols_hu, "b7=b11")$p.value[1],
                     # Framing effects for compound treatment
                     hypotheses(ols_hu, "b4=b8")$p.value[1], 
                     hypotheses(ols_hu, "b4=b12") $p.value[1],
                     hypotheses(ols_hu, "b8=b12")$p.value[1]))

# Switzerland
ch <- data.frame(Switzerland = 
                   # Framing effects for jobs treatment
                   c(hypotheses(ols_ch, "b2=b6")$p.value[1], 
                     hypotheses(ols_ch, "b2=b10")$p.value[1], 
                     hypotheses(ols_ch, "b6=b10")$p.value[1],
                     # Framing effects for prices treatment
                     hypotheses(ols_ch, "b3=b7")$p.value[1], 
                     hypotheses(ols_ch, "b3=b11")$p.value[1], 
                     hypotheses(ols_ch, "b7=b11")$p.value[1],
                     # Framing effects for compound treatment
                     hypotheses(ols_ch, "b4=b8")$p.value[1], 
                     hypotheses(ols_ch, "b4=b12") $p.value[1],
                     hypotheses(ols_ch, "b8=b12")$p.value[1]))

# UK
uk <- data.frame(UK = 
                   # Framing effects for jobs treatment
                   c(hypotheses(ols_uk, "b2=b6")$p.value[1], 
                     hypotheses(ols_uk, "b2=b10")$p.value[1], 
                     hypotheses(ols_uk, "b6=b10")$p.value[1],
                     # Framing effects for prices treatment
                     hypotheses(ols_uk, "b3=b7")$p.value[1], 
                     hypotheses(ols_uk, "b3=b11")$p.value[1], 
                     hypotheses(ols_uk, "b7=b11")$p.value[1],
                     # Framing effects for compound treatment
                     hypotheses(ols_uk, "b4=b8")$p.value[1], 
                     hypotheses(ols_uk, "b4=b12") $p.value[1],
                     hypotheses(ols_uk, "b8=b12")$p.value[1]))

# Table A.2
table <- cbind(de,hu,ch,uk)

p <- xtable(table, label="app:order_support", caption="Test statistic p-values for difference-in-means tests for policy frames")
print.xtable(p, booktabs=TRUE, caption.placement="top")

# Difference-in-means tests for control group
hypotheses(ols_de, "b5=b9")
hypotheses(ols_ch, "b5=b9")
hypotheses(ols_uk, "b5=b9")

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (E) Model specifications with control variables ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define set of control variables
ctrl.soceco <- paste(c("factor(female)", "factor(age_cat)", "factor(edu_cat)", "factor(income_cat)", "factor(wowork)", "factor(retired)"), collapse="+")
ctrl.pol <- paste(c("factor(pol_interest)", "factor(cc_concern)", "factor(lr_cat)"), collapse="+")

# Vectors to reference countries
ktry <- c("germany","hungary","switzerland","uk")
names <- c("Germany", "Hungary", "Switzerland", "UK")

# Loop produces Figures A.2-A.5 and Tables A.3-A.6
for (i in (1:length(ktry))){

# Regression models: collapsed treatment conditions
m1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country==ktry[i],])
m2 <- lm(paste("cbam_support ~ factor(tr_main)", ctrl.soceco, sep="+"), data=dta[dta$country==ktry[i],])
m3 <- lm(paste("cbam_support ~ factor(tr_main)", ctrl.soceco, ctrl.pol, sep="+"), data=dta[dta$country==ktry[i],])
m4 <- lm(paste("cbam_support ~ factor(tr_main) + factor(region)", ctrl.soceco, ctrl.pol, sep="+"), data=dta[dta$country==ktry[i],])

order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  

models.full <- list(
  "Main" = m1,
  "+ w/ soc.econ controls" = m2,
  "+ w/ pol. controls" = m3,
  "+ w/ region FE" = m4)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", facet=TRUE) +  
  labs(x="Coefficients",
       y="Model specifications",
       title=paste("Treatment Effects:",names[i]),
       subtitle="DV: CBAM support (1-5 scale)") +
  xlim(-0.5,0.25) +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF
ggexport(p1, filename=paste0("./Exp_support_controls_",ktry[i],".pdf"), width=10, height=5)

# Produce results table

order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both",
                "factor(female)1"="Female",
                "factor(age_cat)2"="25-34 years old",
                "factor(age_cat)3"="35-44 years old",
                "factor(age_cat)4"="45-54 years old",
                "factor(age_cat)5"="55-64 years old",
                "factor(age_cat)6"="65 years and older",
                "factor(edu_cat)2"="Vocational qualification",                
                "factor(edu_cat)3"="University degree",
                "factor(income_cat)2"="Median income",
                "factor(income_cat)3"="High income",
                "factor(income_cat)4"="Very high income",
                "factor(wowork)1"="Unemployed",
                "factor(retired)1"="Retired",
                "factor(pol_interest)2"="Low political interest",            
                "factor(pol_interest)3"="Medium political interest",
                "factor(pol_interest)4"="High political interest",
                "factor(cc_concern)2"="Low climate concern",
                "factor(cc_concern)3"="Medium climate concern",
                "factor(cc_concern)4"="High climate concern",
                "factor(lr_cat)2"="Central ideology",                 
                "factor(lr_cat)3"="Right ideology"
                )


tab <- modelsummary(models.full,
             statistic = 'conf.int', 
             conf_level = .95,
             coef_rename = order_full,
             coef_omit=c("region"),
             gof_map=c("nobs"),
             title=c(paste0("Regression results for CBAM support with control variables (",names[i],")")),
             output="tinytable") %>%
             theme_tt("placement", latex_float="h!")

tinytable::save_tt(tab, paste0("./main_cbam_controls_",ktry[i],".tex"), overwrite=TRUE)
}



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#                         END OF FILE
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
